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This paper continues an investigation of chaos and chaotic phase mixing in tri- 
axial generalisations of the Dehnen potential which have been proposed to describe 
realistic ellipticals that have a strong density cusp and manifest significant deviations 
^vq , from axisymmetry. Earlier work is extended in three important ways, namely by ex- 

J> ■ ploring systematically the effects of (1) variable axis ratios, (2) 'graininess' associated, 

OO | e.g., with stars and bound substructures, idealised as friction and white noise, and 

0^ ■ (3) large-scale organised motions within a galaxy and a dense cluster environment, 

each presumed to induce near-random forces idealised as coloured noise with a finite 
autocorrelation time. The effects of varying the axis ratio were studied in detail by 
considering two sequences of models with cusp exponent 7 = 1 and, respectively, axis 
ratios a : b : c = 1.00 : 1.00- A : 0.50 and a : b : c = 1.00 : 1.00- A : 1.00 -2A for vari- 
able A. Three important conclusions are (1) that not all the chaos can be attributed to 
• the presence of the cusp, (2) that significant chaos can persist even for axisymmetric 

Qh' systems, and (3) that the introduction of a supermassive black hole can induce both 

q , moderate increases in the relative number of chaotic orbits and substantial increases 

' in the size of the largest Lyapunov exponent. In the absence of any perturbations, 

! the coarse-grained distribution function associated with an initially localised ensem- 

ble of chaotic orbits evolves exponentially towards a nearly time-independent form at 
a rate A that correlates with the typical values of the finite time Lyapunov exponents 
X associated with the evolving orbits. Allowing for discreteness effects and/or an ex- 
ternal environment accelerates phase space transport both by increasing the rate at 
which orbits spread out within a given phase space region and by facilitating diffusion 
along the Arnold web that connects different phase space regions, so as to facilitate an 
approach towards a true equilibrium. The details of the perturbation appear unimpor- 
tant. All that really matters are the amplitude and, for the case of coloured noise, the 
autocorrelation time, i.e., the characteristic time over which the perturbation varies. 
Overall the effects of the perturbations scale logarithmically in both amplitude and 
autocorrelation time. Even comparatively weak perturbations can increase A by a fac- 
tor of three or more, a fact that has potentially significant implications for violent 
relaxation. 
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1 MOTIVATION 

The first paper in this series (Siopis & Kandrup 2000), here- 
after denoted Paper I, began an investigation of phase space 
transport and chaotic phase mixing in triaxial generalisa- 
tions of the Dehnen (1993) potentials which have been pro- 
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posed (e.g., Merritt & Fridman 1996) to model realistic ellip- 
tical galaxies that have a strong density cusp and manifest 
significant deviations from axisymmetry. These correspond 
to potentials associated self-consistently with the mass den- 
sity 

PH = tok m_7(1 + mr(4 " 7) ' (1) 
where 
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As discussed in Paper I, this potential was chosen pri- 
marily to facilitate comparisons with, and to improve the 
relevance of, the results derived here with those of other 
workers who have been using it. It is thought to consti- 
tute a quasi-realistic approximation for the central regions 
of many elliptical galaxies, at least compared with Stackel- 
type potentials which lack a central density cusp. Although 
there is no guarantee that there exist exact self-consistent 
equilibria corresponding to triaxial Dehnen mass distribu- 
tions, insights gained from a study of motion in this po- 
tential can provide clues towards an understanding of more 
generic cuspy triaxial potentials and, importantly, of their 
evolution towards a (quasi-)equilibrium state. 

Paper I explored the effects of varying 7, which controls 
the steepness of the density cusp, and introducing a central 
supermassive black hole of variable mass, as well as allow- 
ing for low amplitude perturbations intended to mimic the 
effects of discreteness effects (i.e., gravitational Rutherford 
scattering) and/or an external environment. However, that 
paper was incomplete in that it did not allow for the effects 
of variable axis ratios, attention focusing exclusively on the 
values c/a = 1/2 and (a 2 - b 2 )/(a 2 - c 2 ) = 1/2 considered 
originally by Merritt & Fridman (1996). 

That work was also incomplete in that the modeling 
of external perturbations was very simplistic. Most of that 
modeling focused on the effects of periodic driving, in which 
the potential is subjected to a time- dependent periodic per- 
turbation characterised by at most three different frequen- 
cies. Although this might seem reasonable when considering 
the effects of a single large companion or a few smaller satel- 
lite galaxies, this is certainly not appropriate to model a rich 
cluster environment, where galaxies tend to be much closer. 
When considering a dense environment one must allow for 
more complex perturbations which, presumably, are far from 
periodic, perturbations which, as described in Paper I, would 
seem better modeled as random kicks of finite duration, i.e., 
coloured noise. Paper I did indeed describe a small num- 
ber of experiments involving coloured noise, but these were 
far from exhaustive. In particular, no effort was made to 
determine the extent to which the detailed form of the per- 
turbation actually matters: all the simulations involved per- 
turbations idealised as an Ornstein-Uhlenbeck process (cf. 
van Kampen 1981). 

The work on chaotic phase mixing described in Paper 
I was also incomplete in the sense that the discussion was 
largely qualitative. It was observed that, as for other poten- 
tials (cf. Kandrup & Mahon 1994, Kandrup 1998b), ensem- 
bles of chaotic orbits typically exhibit a rapid evolution to- 
wards a near-equilibrium (a near- invariant distribution), but 
it was not confirmed explicitly that this evolution proceeds 
exponentially in time. Moreover, there was no systematic 
exploration as to how, for unperturbed systems, the rate A 
at which this evolution proceeds correlates with the size of 
a typical Lyapunov exponent x; ° r , for perturbed systems, 
how A correlates with the amplitude of the perturbation. 

The aim of this second paper is to fill these remaining 
lacunae. Section 2 summarises an investigation of how, ne- 
glecting discreteness and environmental effects but allowing 
for a central supermassive black hole, changes in the axis 
ratios affect (i) the relative number of regular and chaotic 



orbits, (ii) the typical sizes of the Lyapunov exponents, and 
(iii) the overall efficacy of phase space transport. It was 
found that, in general, larger deviations from axisymme- 
try or spherical symmetry tend to increase the fraction of 
chaotic orbits, although considerable chaos can arise even in 
axisymmetric systems, especially at higher energies. Alter- 
natively, strongly triaxial systems do not in general tend to 
have Lyapunov exponents that are much larger than mod- 
erately triaxial systems. The introduction of a black hole 
increases significantly the size of the largest Lyapunov ex- 
ponent, both in absolute units and units of the dynamical 
time to , but does not in general result in a very much larger 
measure of chaotic orbits. The work described here does not 
address the important issue of the relative abundances of dif- 
ferent types of regular orbits. However, this issue is currently 
under investigation in a more general setting (Kandrup & 
Siopis 2003). 

Section 3 focuses on chaotic phase mixing. An analysis 
of coarse-grained distribution functions is used to confirm 
that ensembles of orbits typically evolve exponentially to- 
wards a near-invariant distribution and that, in the absence 
of all perturbations, the rate at which this evolution pro- 
ceeds correlates with the typical size of the finite time Lya- 
punov exponents associated with the ensemble. Discreteness 
effects, modeled (cf. Chandrasckhar 1943) as friction and 
white noise, can dramatically accelerate chaotic phase mix- 
ing both by increasing the rate at which orbits spread out 
within a given chaotic phase space region and by facilitating 
diffusion along the Arnold web that connects different phase 
space regions. 

Section 3 also focuses on how chaotic phase mixing is 
impacted by large-scale organised motions within a galaxy 
or a dense cluster environment, each modeled as coloured 
noise. One principal conclusion of the analysis, consistent 
with Paper I, is that the details of the perturbations seem 
to be largely immaterial: The only things that really seem 
to matter are (1) the amplitude and (2) the autocorrelation 
time (i.e., the characteristic time scale associated with in- 
dividual kicks), both of which can be readily estimated via 
dimensional analysis. In other words, the effects of the envi- 
ronment seem to be insensitive to details which are difficult 
to ascertain observationally! The other principal conclusion 
is that choices of amplitude and autocorrelation time appro- 
priate for real galaxies can in fact lead to significant effects 
on time scales short compared with the age of the Universe. 

Section 4 summarises the principal conclusions and Sec- 
tion 5 discusses potential implications for real galaxies. 

2 CHAOS AS A FUNCTION OF AXIS RATIO 
AND BLACK HOLE MASS 

2.1 The numerical experiments 

In order to understand the effects of variable axis ratio, two 
different sequences of models were considered in detail. The 
first involved sweeping through a variety of triaxial con- 
figurations which connected prolate and oblate spheroids. 
Specifically, this entailed assuming axis ratios a : b : c — 
1.00 : 1.00 — A : 0.50, with A allowed to vary between 
0.00 and 0.50 by increments S — 0.05, but also consider- 
ing several extra models with A closer to 0.00 and 0.50. 
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The second sequence involved triaxial deviations from a 
spherical system, corresponding to axis ratios a : b : c = 
1.00 : 1.00 - A : 1.00 - 2A. 

Attention was restricted primarily to models with cusp 
index 7=1 but three different black hole masses were con- 
sidered, namely M B h = 0, M B h = 10~ 3 , and M B h = 10~ 2 . 

The principal focus was on determining the statisti- 
cal properties of orbit ensembles as a function of energy 
E. This was done by considering ten different energies, 
— 1.0 < E < —0.1, and, for each energy, selecting 1000 'rep- 
resentative' initial conditions. (The model with axis ratio 
1.00 : 0.75 : 0.50, which arises in both sequences, was stud- 
ied for two different sets of 1000 initial conditions and it 
was confirmed that the conclusions were in agreement sta- 
tistically.) Arguably the most honest sampling of any given 
energy entails a uniform sampling of the constant energy 
hypersurface; and, for this reason, orbits were selected to 
sample the microcanonical distribution 

H<x6 D {H -E), (3) 

with H the Hamiltonian. This distribution turns out to be 
difficult and expensive to sample directly, especially since 
the potential V cannot be expressed analytically. For this 
reason an indirect approach was used. By integrating over 
the dependence on velocity, it is easily seen that a micro- 
canonical distribution corresponds to a configuration space 
distribution 

f (E - V) 1 ' 2 if V(r) < E; 
/(r) oc \ (4) 
I if V(r) > E. 

Alternatively, the velocity distribution at any given config- 
uratio n space is isotropic, with each mass having a speed 
v — \fl{E — V). To obtain a random, uniform sampling of 
the constant energy surface, it therefore suffices to (1) sam- 
ple /(r) to generate a collection of 1000 configuration space 
points and (2) to each point assign a velocity of magnitude 
v — \[%E — V) oriented in a randomly chosen direction. 

Each orbit was integrated for a time > 200to using a 
variable timestep integrator with accuracy parameter 10~ 8 
which, in every case, conserved energy to at least 1 part in 
10 5 . Estimates of the largest (finite time) Lyapunov expo- 
nent were obtained by tracking the evolution of a nearby 
orbit which was periodically renormalised in the usual way. 
The dynamical time to for a given choice of energy and axis 
ratios was identified as follows: For the orbits in each ensem- 
ble, define t x - and, analogously, the quantities t y and t z - 
as the mean time between successive crossings of the x — 
plane; and, given t x , t y , and t z , define 

t D = 2n(t x + t y +t z ). (5) 

This definition seems well motivated physically; the choice 
of normalisation factor 2n ensured that, for the 'maximally 
triaxial' model first considered in detail by Merritt and Frid- 
man (1996), this definition of to agreed with the Merritt- 
Fridman definition for all energies to within 3%. 

The degree of chaos manifested by the different ensem- 
bles was quantified using two complementary diagnostics, 
namely finite time Lyapunov exponents (Grassberger et al 
1988), which probe the degree of exponential sensitivity ex- 
hibited by different orbits, and orbital complexities (Kan- 
drup, Eckstein, and Bradley 1997), which probe the extent 



to which the power associated with an individual orbit is 
concentrated at or near a few special frequencies. This en- 
tailed determining for each orbit the quantities n x , n y , and 
n z , defined, respectively, as the minimum number of fre- 
quencies required to capture a fixed fraction k of the power 
in each direction, and then assigning a total complexity 

n = n x + n y + n z . (6) 

Given temporal discreteness effects reflecting the fact that 
the orbital data were sampled at intervals ~ 0.25£d, the 
cleanest distinctions between different orbits were obtained 
for k w 0.9. 

The degree to which individual chaotic orbits are 
'sticky' (Contopoulos 1971), i.e., that they can remain 
'stuck' in a small part of the accessible phase space for a 
comparatively long time, was probed by computing finite 
time Lyapunov exponents for long time integrations of sev- 
eral different initial conditions on the same phase space hy- 
persurface and determining the time scale (or, in some cases, 
a lower bound) on which the exponents for different orbits 
converge towards a single asymptotic value. 

Because it can take an orbit a very long time to uni- 
formly sample the accessible phase space, it can be difficult 
and expensive computationally to compute estimates of the 
largest Lyapunov exponent as a function of axis ratio, black 
hole mass, and energy. For this reason, such estimates were 
instead obtained typically by computing the mean value of 
the finite time Lyapunov exponents for the chaotic orbits 
in a given 1000 orbit ensemble, a procedure which is jus- 
tified theoretically given the assumption of ergodicity. For 
the case of models where only a very few orbits are chaotic, 
this procedure is suspect statistically, so that the value of 
the largest Lyapunov was computed by averaging over the 
results of extremely long time integrations of four chaotic 
initial conditions. 

Uncertainties in the relative measure of chaotic or- 
bits were estimated assuming N 1 ^ 2 statistics. Uncertain- 
ties in the estimated Lyapunov exponent were obtained by 
analysing separately the two halves of the 1000 orbit ensem- 
bles. 



2.2 Results 

Consider first the sequence of models extending between 
oblate and prolate axisymmetric configurations. Here the 
most obvious point is that, overall, triaxial models which 
manifest larger deviations from axisymmetry tend to admit 
larger measures of chaotic orbits. This is, e.g. evident from 
FIGURE 1, which exhibits the fraction / of chaotic orbits 
as a function of the intermediate axis b for ten different en- 
ergies. 

Nevertheless, as is evident from FIGURE 1, even ax- 
isymmetric configurations can admit significant measures of 
chaos, at least at higher energies. This chaos appears to re- 
flect the large scale structure of the bulk potential, not the 
presence of a central cusp. Indeed, for very high energies, 
one finds virtually identical numbers for a cuspless 7 = 
model although, for lower energies, the 7 = admits a much 
smaller measure of chaotic orbits. This is illustrated in FIG- 
URE 2 which, for a model with axis ratio 1.00 : 0.50 : 0.50, 
exhibits the ordered values of the 1000 computed Lyapunov 
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Figure 1. (a) The fraction of orbits with energy E = —1.0 that 
are chaotic for models with axis ratios 1.00 : b : 0.50 and black 
hole mass Mbh = (solid curve), Mbh = 10 -3 (dashed curve), 
and Mbh = 10~ 2 (dot-dashed curve) (b) The same for E = —0.9. 
(c) The same for E = -0.8. (d) The same for E = -0.7. (c) The 
same for E = —0.6. (f) The same for E = —0.5. (g) The same 
for E = -0.4. (h) The same for E = -0.3. (i) The same for 
E = -0.2. (j) The same for E = -0.1. 



exponents for both 7 = 1 and 7 = for energies ranging be- 
tween E — —0.1 and E = —0.5. However, despite this chaos 
at higher energies, except for the case of prolate axisym- 
metric configurations there tends overall to be more chaos 
at lower energies, presumably associated with the cusp: this 
chaos is, e.g., absent for analogous models with 7 = 0. For 
a broad range of parameter values, the measure of chaotic 
orbits ranges between about 10% and 50%. 

That chaos was observed in the prolate model with 
axis ratio 1.00 : 0.50 : 0.50 but not in the oblate model 
with 1.00 : 1.00 : 0.50 raises the question as to whether 
chaos is in fact generic in the axisymmetric Dehnen po- 
tentials. This issue was addressed by considering a variety 
of prolate and oblate models with, respectively, axis ratios 
1.00 : 1.00 - A : 1.00 - A and 1.00 : 1.00 : 1.00 - A. The 
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Figure 2. (a) An ordered plot of the computed Lyapunov expo- 
nents for the 7 = 1 model with a : b : c = 1.00 : 0.50 : 0.50 and 
energy E = —0.1. (b) The same for 7 = 0. (c) and (d) The same 
for E = -0.2. (c) and (f) The same for E = -0.3. (g) and (h) 
The same for E = —0.4. (i) and (j) The same for E = —0.5. 



net result is that chaos can arise for both prolate and oblate 
systems, but that, for the case of oblate systems, a substan- 
tially larger deviation from sphericality is required. For the 
sequence with 1.00 : 1.00 — A : 1.00 — A, a significant amount 
of chaos, ~ 2%, is observed already at energy E = —0.1 for 
A = 0.2; for the sequence with 1.00 : 1.00 : 1.00 - A, one 
requires a value as large as A = 0.6 to find a comparable 
amount of chaos. This is illustrated in FIGURE 3, the top 
two panels of which exhibit the fraction / of chaotic orbits 
with energy E = —0.1 for different choices of axis ratio. 
The lower two panels exhibit estimates of the largest Lya- 
punov exponent for the same configurations. In each case, 
the onset of chaos is characterised by very small Lyapunov 
exponents which, however, increases monotonically or near- 
monotonically with increasing A. 

It is interesting that an analogous result has been ob- 
tained for the so-called thermal equilibrium model, a stan- 
dard pedagogical example (cf. Brown & Reiser 1995) from 



© 2002 RAS, MNRAS 000, 1-?? 



Chaos and chaotic phase mixing in cuspy triaxial potentials 5 




0.2 




4 0.2 



0.052 — 
0.039 - 
0.026 - 

0.01 3 - 
0.000 L 
1 .0 



(c) 




Am*/ 

J 5 



0.8 0.6 0.4 0.2 



0.052 
0.039 

0.013 ■ 
0.000 : 
1 .0 



(d) 



0.8 0.6 0.4 0.2 



Figure 3. (a) The fraction / of chaotic orbits for E = —0.1 and 
variable axis ratios a : b : c = 1.00 : c : c. (b) The same for 
a : b : c = 1.00 : 1.00 : c. (c) The largest Lyapunov exponent for 
the parameter values in (a), (d) The same (b). 



the physics of charged particle beams. For this system, which 
corresponds to a self-interacting nonneutral plasma in ther- 
mal equilibrium confined by an anisotropic harmonic oscil- 
lator potential, one finds (Bohn & Sideris 2003) that even 
slightly prolate configurations tend to admit large amounts 
of chaos, whereas comparably aspherical oblate configura- 
tions exhibit little, if any, chaos. 

The remaining interesting point is that, even when 
an axisymmetric model admits no chaotic orbits, relatively 
small perturbations away from axisymmetry suffice to trig- 
ger a significant amount of chaos, 10% or more. This is in 
qualitative agreement with recent computations by El-Zant 
and Shlosman (2002), who explored the effects of bars with 
variable amplitude on orbits in an otherwise axisymmetric 
potential. 

Overall, when scaled in physical units, the size of the 
largest Lyapunov exponent xto ~ 0.2, i.e., the growth time 
is roughly 5trj. This is, e.g. evident from FIGURE 4 which 
exhibits estimates of the largest Lyapunov exponent for the 
same ensembles used to generate FIGURE 1. For lower ener- 
gies, this conclusion is again comparatively insensitive to the 
choice of axis ratio. However, the behaviour at much higher 
energies, where chaos is triggered by the bulk potential, is 
different. In this case, the value of x^d is actually maximised 
for the axisymmetric model with TOO : 0.50 : 0.50 and de- 
creases monotonically as one moves along the sequence to- 
wards TOO : TOO : 0.50. 

The introduction of a black hole has only a minimal 
effect on high energy orbits which spend little time in the 
central region where they can 'feel' the gravitational influ- 
ence of the hole. At lower energies, however, the presence 
of a black hole tends to increase the relative measure of 
chaotic orbits, albeit not by all that much. Thus, e.g. as is 
evident from FIGURE 1, the introduction of a hole with 
Mbh = 10~ 2 invariably increases the fraction / by 50% or 
less. There always remains a significant measure of regular 
orbits. Interesting also is the fact that curves of /(A) for 
models with and without a black hole manifest similar cur- 
vatures: the only obvious difference is that, for Mbh = 0, 
there are somewhat fewer chaotic orbits. 

Although a black hole changes / by only a relatively 
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Figure 4. (a) Estimates of the largest Lyapunov exponent for 
chaotic orbits with energy E = —1.0 for models with axis ratio 
1.00 : b : 0.50 and black hole mass Mbh = (solid curve), 
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curve) (b) The same for E = —0.9. (c) The same for E = —0.8. 
(d) The same for E = -0.7. (e) The same for E = -0.6. (f) The 
same for E = —0.5. (g) The same for E = —0.4. (h) The same for 
E = -0.3. (i) The same for E = -0.2. (j) The same for E = -0.1. 



small amount, it typically occasions a substantial increase in 
the value of the largest Lyapunov exponent for low energies. 
Moreover, one sees that although / appears to vary smoothly 
with shape, the value of the largest Lyapunov exponent can 
exhibit a more complex dependence on A. To summarise: 
the introduction of a black hole typically does not result in 
a huge increase in the number of chaotic orbits, but it does 
increase significantly the degree of exponential sensitivity 
exhibited by those orbits which are chaotic. 

A consideration of the sequence which starts from 
spherical and becomes triaxial also yields several interest- 
ing conclusions. Most obvious from FIGURE 5 is the fact 
that, for all energies, the relative measure of chaotic orbits 
is a monotonically increasing function of A, i.e., the devia- 
tion from sphericality. Even a deviation as small as A = 0.01 
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Figure 5. (a) The fraction of orbits with energy E = —1.0 that 
are chaotic for models with axis ratios 1.00 : b : c for b = 1 — A 
and c = 1 — 2A, and black hole mass Mbh = (solid curve), 
Mbh = 10~ 3 (dashed curve), and Mbh = 10 -2 (dot-dashed 
curve) (b) The same for E = —0.9. (c) The same for E = —0.8. 
(d) The same for E = -0.7. (c) The same for E = -0.6. (f) 
The same for E = -0.5. (g) The same for E = -0.4. (h) The 
same for E = —0.3. (i) The same for E = —0.2. (j) The same 
for E = —0.1. The unusual appearance of panel (a) reflects the 
fact that, for nearly spherical models, the minimum energy is 
somewhat larger than E = —1.0. 



suffices to trigger a nonzero measure of chaotic orbits. How- 
ever, a comparatively large A is required to trigger as much 
as (say) 20% chaotic orbits. 

Overall, as is evident from FIGURE 6, at least for 
Mbh = the value of the largest Lyapunov exponent also 
appears to be a monotonically increasing function of A. 
However, especially for lower energies the increase is rela- 
tively minimal for A > 0.1 or so. In other words, larger de- 
viations from sphericality yield a larger measure of chaotic 
orbits, but the degree of exponential sensitivity, as probed 
by the largest Lyapunov exponent, does not change all that 
much. 
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Figure 6. (a) Estimates of the largest Lyapunov exponent for 
chaotic orbits with energy E = —1.0 for models with axis ratio 
1.00 : b : c for b = 1 - A, and c = 1 - 2A, and black hole 



mass Mbh = (solid curve), Mbh 



10 



-3 



(dashed curve), and 



M B h = 10~ 2 (dot-dashed curve) (b) The same for E = -0.9. 
(c) The same for E = -0.8. (d) The same for E = -0.7. (c) 
The same for E = -0.6. (f) The same for E = -0.5. (g) The 
same for E = —0.4. (h) The same for E = —0.3. (i) The same for 
E = -0.2. (j) The same for E = -0.1. 



As for the other sequence, the addition of a black hole 
again tends to increase the abundance of chaotic orbits, al- 
beit not by all that much. Adding a black hole also leads 
to larger Lyapunov exponents but, unlike the case when 
Mbh = 0, this exponent is not a strictly monotonically 
increasing function of A. For the larger black hole mass, 
Mbh = 10~ 2 , the largest exponent increases rapidly with 
A up to a value A ~ 0.1 — 0.2 but then becomes a compar- 
atively flat function of A. 

For virtually all choices of axis ratio, some chaotic orbit 
segments were so nearly regular that it was highly nontrivial 
to distinguish then from regular segments. (For this reason 
distinctions between regular and chaotic were made by us- 
ing both Lyapunov exponents and complexities as comple- 
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Figure 7. (a) Time-dependent estimates of the true (cumula- 
tive) Lyapunov exponent for two orbits with i? = —1.0 evolved in 
the absence of a black hole in a model with axis ratio 1.00 : 0.95 : 
0.50. (b) The same for E = —0.1. (c) and (e). Estimates of finite 
time Lyapunov exponents for the orbits exhibited in (a), gener- 
ated by partitioning the data into a large numer of segments and 
analysing them individually, (d) and (f) The same for the orbits 
in (b). 



mentary diagnostics.) Moreover, the structure of the chaotic 
phase space regions is often complex in the sense that two 
chaotic segments evolved in the same potential with the 
same energy can have finite time Lyapunov exponents with 
significantly different values for surprisingly long times. 

This sort of stickiness appears to be especially pro- 
nounced for systems that are comparatively close to axisym- 
metric. One example is exhibited in FIGURE 7, which was 
generated from orbit segments evolved in a model with axis 
ratio 1.00 : 0.95 : 0.50. Each of the top two panels shows 
the time-dependent x{t)i generated in the usual way by also 
tracking the evolution of a nearby orbit periodically renor- 
malised, for two orbits each with energies E — —1.0 and 
E — —0.1 evolved for a total time t — 32768ir>. In each case, 
estimates of the largest Lyapunov exponent were computed 
at intervals of lto- The fact that the curves are not decay- 
ing uniformly towards zero is prima facia evidence that the 
orbits are not regular. How the actual degree of exponential 
sensitivity varies with time can be gauged from the lower 
four panels, which exhibit finite time Lyapunov exponents 
for segments of these orbits. These panels were constructed 
by extracting from each orbit a collection of 32768 finite time 
Lyapunov exponents {x(li)} f° r intervals U < t < U + 1 
(i = 0, 32767) and then smoothing the resulting data by 
performing a boxcar average over 400 adjacent points. It is 
evident that, in each case, these finite time exponents exhibit 
considerable variability. However, it is also clear that the or- 
bits for which the cumultative x{t) is smaller tend to spend 
a considerable amount of time in phase space regions where 
the degree of exponential sensitivity is very small, whereas 



the orbits with larger \ tend systematically to avoid these 
regions. 



3 CHAOTIC PHASE MIXING 
3.1 The numerical experiments 

As in Paper I, the experiments involving chaotic phase mix- 
ing entailed (i) selecting localised ensembles of 1600 initial 
conditions all with the same energy E, (ii) evolving each 
orbit in the ensemble for a time t = 200tc, (iii) determin- 
ing the extent to which the evolving ensemble exhibited a 
coarse-grained evolution towards a nearly time-independent 
state, i.e., a near-invariant distribution, and (iv) determin- 
ing how this evolution was affected by discreteness effects 
and/or an external environment, modeled as friction and 
white or coloured noise. Attention focused primarily on weak 
perturbations, characterised typically by a relaxation time 
tR much longer than the time scales of interest, so that, in 
a first approximation, the energies of individual orbits are 
nearly conserved and one can view the perturbations as sim- 
ply accelerating diffusion on constant energy hypersurfaces 
(compare Weinberg 2001a,b). 

A real galaxy is of course characterised by a complex, 
many-body potential; and, for that reason, it is not com- 
pletely obvious that constructs like cantori (Mather 1982) 
or the Arnold web (Arnold 1964), which can be proven rig- 
orously for smooth potentials, are applicable to galaxies (see, 
e.g., Contopoulos 2002 for a pedagogical discussion of can- 
tori and the Arnold web). However, recent work (Sideris & 
Kandrup 2002) suggests strongly that, at least for the case of 
systems in a time-independent equilibrium, such constructs 
do remain applicable. Specifically, a comparison of the same 
initial conditions evolved both in the smooth potential and 
in fixed (in space and time) iV-body realisations of the cor- 
responding density distribution reveals that, both for or- 
bit ensembles and for individual orbits, iV-body trajectories 
are extremely well modeled by smooth potential orbits per- 
turbed by Gaussian white noise with amplitude consistent 
with that predicted by a Fokker-Planck description (Chan- 
drasekhar 1943). Indeed, it is possible to extract estimates 
of finite time Lyapunov exponents for the smooth potential 
from such iV-body trajectories (Kandrup & Sideris 2003). 

This suggests that, in a first approximation, one can 
interpret orbits in JV-body systems as exhibiting (Lichten- 
berg & Lieberman 1992) 'intrinsic diffusion', just as in the 
smooth potential, so that, e.g., one can visualise orbits pass- 
ing through cantori in the fashion described by the 'turnstile 
model' of MacKay, Meiss, and Percival (1984). In a next ap- 
proximation, it would then seem reasonable to model dis- 
creteness effects using Fokker-Planck or Langevin simula- 
tions (Chandrasekhar 1943), the friction and noise associ- 
ated with such simulations being interpreted (Lichtenberg 
& Lieberman 1992) as a source of 'extrinsic diffusion' which 
tends, generically, to accelerate phase space transport. 



3.1.1 White noise 

The initial conditions were generated by sampling a tiny 
phase space hypercube of characteristic size r ~ 0.01 or less. 
In the absence of perturbations, orbits were generated from 
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these initial conditions by solving the Hamilton equations 
appropriate for motion in the potential (1), using a vari- 
able time step integrator which typically conserved energy 
to better than one part in 10 5 . Following Chandrasekhar 
(1943), discreteness effects, i.e., the effects of gravitational 
Rutherford scattering, were modeled as resulting in fric- 
tion and additive white noise connected by a Fluctuation- 
Dissipation Theorem. The evolution of individual orbits thus 
entailed solving a Langevin equation (Chandrasekhar 1943, 
van Kampen 1981) of the form 



dV(r) 



-nv a + F a , (a = x,y,z), 



(7) 



dt 2 dx a 

with n a constant coefficient of dynamical friction and F 
a 'stochastic force.' Assuming in the usual fashion that F 
corresponds to homogeneous Gaussian noise, its statistical 
properties are characterised completely by its first two mo- 
ments, which take the form 



(F„(t))=0 and 

{FaitJFifo)) = 6 ab K(t! 



■t 2 ), (a,b = x,y,z). 



(8) 



The assumption that the noise be white implies further that 
the autocorrelation function K is proportional to a Dirac 
delta, so that 



K(t) = 2t)&5 d (t). 



(9) 



Physical considerations dictate that the 'temperature' 
~ \E\, so the simulations were all performed assuming 
O = —E. In this case, n defines the relaxation time tn 
on which the energy of an orbit evolved in an otherwise 
fixed potential would change significantly: tn = r)^ 1 The 
Langevin equation was solved using an algorithm developed 
by Griner, Strittmatter, & Honerkamp (1988) (see also Hon- 
erkamp 1994). 

Orbital data, recorded at intervals St = O.lto were 
binned into rectangular grids comprised of 20 x 20 cells 
so as to construct coarse-grained distribution functions 
f{Z a ,Zb,i) for pairs of phase space variables, i.e., aj^b = 
x, y, z, v x , v y , v z . An examination of successive 'snapshots' 
revealed a systematic tendency for the ensemble to disperse 
and, eventually, to approach a nearly time-independent 
state. For this reason, the last 500 snapshots, appropriate 
for 150. lto < t < 200. 0to were combined to generate a nu- 
merical representation of a coarse-grained near-invariant dis- 
tribution, i.e., 



ZUUU 

- wo E '<*>■ 



(10) 



The approach of f(t) towards the near- invariant f„i v was 
quantified by computing an L 2 "distance" between f(t) and 
f niv via the natural prescription (cf. Kandrup & Mahon 
1994, Merritt & Valluri 1996, Kandrup 1998b) 

/ \ 1/2 

Z)q Z)b \f(Za,Zb,t)-f niv (Zg,Z b )[ 



Df(Z a ,Z b ,t) 



(11) 



3.1.2 Coloured noise 

Discreteness effects such as stellar encounters correspond to 
instantaneous random kicks, and can be adequately mod- 
eled as white noise. However, other discreteness effects (e.g., 



large-scale organised motions inside a galaxy or perturba- 
tions caused by an external environment) are better mod- 
eled as random kicks of finite dutation. Assuming that these 
kicks constitute a random Gaussian process, one is led once 
again to solve the Langevin equation (7), the only differ- 
ence being that the autocorrelation function K(ti — £2) is 
no longer idealised as a Dirac delta ('coloured noise'). 

Two specific choices for the autocorrelation function 
were considered. One corresponds to a so-called Ornstein- 
Uhlenbeck process (cf. van Kampen 1981), for which K de- 
creases exponentially in time, i.e., 



K(t) = anO exp(— a\r\). 



(12) 



This process is characterised by three parameters: r\ and 
O have the same meaning as for the case of white noise, 
whereas the autocorrelation time 



!~TK{T)dT 

J-K(r)dr 



(13) 



sets the time scale on which the random forces change ap- 
preciably. The normalisation in eq. (12) ensures that the 
diffusion constant D that would enter into a Fokker-Planck 
description, 



drK{r) = 2®n, 



(14) 



is independent of a. 

The other choice (cf. Pogorelov & Kandrup 1999) cor- 
responds to an autocorrelation function of the form 

K{r) = ^ exp(-a|r|) (l + a\r\ + ^r 2 ) . (15) 

For fixed a this autocorrelation function decays somewhat 
more slowly, the autocorrelation time now equaling t c = 
2/q. In both cases, the integrations were performed using a 
variant of an algorithm summarised in Pogorelov & Kandrup 
(1999). White noise can be viewed as a singular a —> 00 limit 
of either of these two processes. 



3.2 Unperturbed Hamiltonian evolution 

As for simpler potentials (Kandrup 1998a), an initially lo- 
calised ensemble of chaotic orbits tends to disperse expo- 
nentially at a rate A comparable to the typical size of the 
finite time Lyapunov exponents, x> for the ensemble. This 
implies, e.g., that quantities like a x and o~ Vx , the disper- 
sions in position and velocity associated with the ensemble, 
initially grow exponentially. More significant, however, for 
the problem of chaotic phase mixing, is the fact that this 
evolution also entails a comparatively efficient approach to- 
wards a near- invariant distribution f n i v . In most cases, this 
evolution towards f„i V is approximately exponential, so that 
Df(t) oc exp(— At); and, in agreement with Kandrup & Ma- 
hon (1994) and Merritt & Valluri (1996), one finds typically 
that A is again comparable in magnitude to, albeit some- 
what smaller than, a typical finite time Lyapunov exponent 
X- 

A typical example of this behaviour is illustrated in 
the top three left hand panels of FIGURE 8, which exhibit 
Df(x,y), Df(y,z), and Df(z,x) for one initially localised 
ensemble of chaotic orbits in the lowest energy shell evolved 
in the triaxial Dehnen potential with axis ratios c/a — 1/2 
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Figure 8. (a) The L distance Df(x,y,t) between f(x,y,t) and 
a near-invariant f n iv(%,y) computed for an initially localised 
ensemble of 1600 chaotic orbits in the lowest shell evolved in 
the triaxial Dchncn potential with axis ratios c/a = 1/2 and 
(a 2 — b 2 )/(a? — c 2 ) = 1/2. (b) Df(x,y,t) for another ensemble 
evolved in the same potential with the same energy, (c) Df(y, z, t) 
for the ensemble in (a), (d) Df(y, z, t) for the ensemble in (b). (e) 
Df(z, x, i) for the ensemble in (a), (f) Df(z, x, t) for the ensemble 
in (b). (g) JV[x], the distribution of finite time Lyapunov expo- 
nents for the initial conditions in (a) evolved for a time t = 200t o ■ 
(h) N[x] for the initial conditions in (b) evolved identically. 



and (a 2 - 6 2 )/(a 2 - c 2 ) = 1/2. The bottom left hand panel 
exhibits a distribution of finite time Lyapunov exponents for 
the same ensemble computed for a total interval t = 200izj. 
It should be noted that the saturation of In Df at a value 
~ — 3.5 observed in these panels is a numerical artifact, 
rather than a real physical effect: Even if the data points 
used to generate f(t) and /„»„ had been obtained by ran- 
domly sampling exactly the same continous distribution, the 
computed distance Df between them would be nonzero be- 
cause of finite number statistics. 

Occasionally, however, one finds that the approach to- 
wards f n iv is significantly less efficient. This is, e.g., illus- 
trated by the right hand panels of FIGURE 8, which exhibit 
Df and the distribution of finite time Lyapunov exponents 
for a different ensemble of chaotic orbits evolved in the same 
potential and with the same energy. At very early times, the 
distances Df computed for this ensemble are not apprecia- 
bly different from those observed for the first ensemble, but 
after t ~ 15ir> or so, it is clear that the approach towards a 
near- invariant /„»„ has become appreciably less efficient. 
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Figure 9. (a) The dispersion cr x for the ensemble exhibited in the 
left panels of FIG. 1. (b) u x for the ensemble in the right panels, 
(c) and (d) The corresponding a y 's. (c) and (f) The corresponding 
a z 's. 



As is evident from FIGURE 9, the qualitative differ- 
ences exhibited between these two ensembles are also re- 
flected by the evolution of quantities like the dispersions a x , 
a y , and a z . At very early times the dispersions for the two 
different ensembles evolve in a comparatively similar fash- 
ion, but three obvious differences are evident at later times: 
The later time evolution for the right hand ensemble is sig- 
nificantly less smooth, the dispersions exhibiting high fre- 
quency variability which only damps over times t > WOto- 
Moreover, it is apparent that, even after a time as long as 
t = 200to, o z is exhibiting a systematic secular evolution. 
Both these features would suggest that the ensemble is com- 
paratively inefficient in achieving a near-invariant distribu- 
tion, consistent with what one infers from FIGURE 8. The 
third point is that the asymptotic values towards which the 
tTj converge are different for the two ensembles. This reflects 
the fact that the near- invariant fniv's for the two ensem- 
bles differ. It would appear that, at early times, the two 
ensembles are restricted to different phase space rehions al- 
though, at late times, they will presumably converge towards 
the same invariant f inv . 



3.3 Evolution including discreteness effects 
modeled as white noise 

Discreteness effects, modeled as friction and white noise, 
can accelerate the approach towards an invariant or near- 
invariant distribution in two relatively distinct fashions. 
On the one hand, such irregularities can accelerate phase 
space transport between different chaotic phase space re- 
gions which, presumably, are connected by cantori or the 
Arnold web. On the other, they can significantly increase the 
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efficiency with which orbits disperse within a single chaotic 
phase space region so that, e.g., the rate A associated with 
the initial approach towards a near-invariant distribution 
becomes larger. 

That discreteness effects can accelerate phase space 
transport through topological obstructions is a fact that has 
been recognised previously, both in the context of the triax- 
ial Dehnen potentials (Siopis & Kandrup 2000) and for sim- 
pler two and three degree of freedom systems (cf. Kandrup, 
Pogorelov & Sideris 1999). As first discussed in Pogorelov 
& Kandrup (1999), this can be understood as reflecting the 
fact that weak irregularities wiggle the orbits, thus assisting 
them in finding phase space holes. 

That discreteness effects can also accelerate phase space 
transport within a given phase space region does not seem to 
have been recognised previously. The fact that introducing 
weak irregularities can increase A is significant both for its 
physical implications - which will be discussed more care- 
fully below - and for what it may suggest about the physical 
mechanism responsible for chaotic phase mixing. It is hardly 
surprising that an initially localised ensemble of chaotic or- 
bits evolved into the future will disperse at a rate A compa- 
rable to a typical finite time Lyapunov exponent \ f° r the 
orbits in the ensemble. However, it is not completely obvi- 
ous how, if at all, the rate A associated with the exponential 
approach towards a near- invariant distribution should cor- 
relate with x- It is thus interesting, and significant, that, 
for a variety of different potentials (cf. Kandrup & Mahon 
1994, Merritt & Valluri 1996), in the absence of perturba- 
tions there seems to be a strong correlation between \ and 
A, although the correlation is not completely linear. Overall, 
A seems to be a factor of 2 or 3 smaller than x- 

Except for the very largest perturbations, r\ > 10~ 2 ' 5 
or so, friction and noise do not affect the exponential rate 
at which orbits in an ensemble disperse, although they do 
have a subexponential effect. If, e.g., one selects an ensemble 
that samples a very tiny region (or, as an extreme case, 
simply tracks multiple noisy realisations of the same initial 
condition), the dispersions will grow as 

a x tx (Or]) 1 ' 2 exp(xi) (16) 

(cf. Habib, Kandrup, & Mahon 1997). The irregularities 
facilitate the spreading of the ensemble, but they do not 
change the overall rate! By contrast, friction and noise do 
increase the rate A at which the orbit ensemble evolves to- 
wards a near-invariant distribution. 

This increase is evident visually from FIGURES 10 and 
11, which exhibit Df(x,y,t) for the same two sets of initial 
conditions used to generate FIGURES 8 and 9, now evolved 
allowing for friction and white noise with variable tR be- 
tween 10 4 ir> and 10 6,5 tn. Best fit values of A for the initial 
interval < t < 6io, both for these and other values of tR, 
are exhibited in the top two panels of FIGURE 12. In each 
case, the computed value averages over the fifteen possible 
choices of f(Z a ,Z b ), and the formal error bars correspond 
to the standard deviations associated with the mean. The 
dashed lines correspond to the best fit value of A for simu- 
lations without any friction or noise. 

Although a useful probe of the early stages of the evolu- 
tion towards a near-invariant distribution, A misses the clear 
differences between the two ensembles that are obvious vi- 
sually in FIGURE 8. An alternative probe, more sensitive 



to somewhat later time evolution, is the time r required be- 
fore Df(t) becomes smaller than some fiducial value. The 
bottom two panels of FIGURE 12 exhibit r.05, the time re- 
quired for Df(t), which starts from an initial value Df = 1, 
to decrease to a value Df(t) — 0.05. Perhaps the most strik- 
ing thing about FIGURE 12 is the fact that even very weak 
friction and noise, corresponding to tR > 10 6 ir>, is enough 
to dramatically accelerate the approach of the second en- 
semble towards a near-invariant distribution. The best fit 
values of A in the absence of noise and for t R = W &5 t D are 
virtually identical, but the noisy simulation is characterised 
by a value of r.05 that is less than half as large as the value 
assumed in the absence of noise! 

Another interesting feature, again evident from FIG- 
URE 12, is that A and T.05 both exhibit a roughly logarith- 
mic dependence on r\ and tR. This logarithmic dependence, 
observed also when probing the effects of friction and noise 
on phase space transport through cantori in two-degree-of- 
freedom systems (Pogorelov & Kandrup 1999), implies that 
the effects of the perturbations only turn on gradually. There 
is no critical threshhold amplitude above which the pertur- 
bations immediately become important. Equation (12) fa- 
cilitates at least a heuristic explanation of why the time r 
required to converge towards f niv scales logarithmically in 
n or tR. If one assumes, simplistically, that r should scale as 
the time required for the ensemble to expand to a size com- 
parable to the accessible phase space region, r corresponds 
to a time when the configuration space dispersions have as- 
sumed some fiducial value. Assuming, however, that this be 
true, eq. (16) implies that r oc const + \11tR. 

The effects of noise on orbits inside a given phase space 
region also resemble the effects of noise on diffusion between 
different phase space regions in one other important respect: 
the details seem largely unimportant. Turning off the fric- 
tion but retaining the noise, or making the white noise mul- 
tiplicative (i.e., state-dependent), so that n is a nontrivial 
function of x and/or v, seems largely irrelevant. 

In any event, what is apparent, e.g., from FIGURE 12, 
is that even comparatively low amplitude perturbations can 
dramatically accelerate chaotic phase mixing within a given 
phase space region. Friction and noise corresponding to a 
relaxation time as long as tR ~ 10 5 tc can increase the rate 
of chaotic phase mixing by ~ 50%; friction corresponding to 
t R ~ 10 4 t D can increase A by a factor of two. What this sug- 
gests is that, even in settings where diffusion through cantori 
is comparatively unimportant, friction and noise can play an 
important role in enhancing the overall efficacy with which 
orbits disperse and, presumably, the rate at which a config- 
uration displaced from equilibrium can readjust towards a 
new equilibrium. 



3.4 Evolution including perturbations modeled as 
coloured noise 

Just as for white noise, it was found that coloured noise can 
significantly accelerate phase space transport both by facil- 
itating transport along the Arnold web and by enhancing 
diffusion within a single chaotic phase space region. In par- 
ticular, one discovers once again that the initial evolution 
towards a near-invariant distribution f„i V is typically well 
fit by an exponential, at least at early times, and, for fixed 
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Figure 10. The L 2 distance Df(x,y,t) between f(x,y,i) and a 
near-invariant fi nv (x,y) computed for the first ensemble in FfG. 
1, now perturbed by friction and white noise with = — E and 
variable r) = t" 1 . (a) t R = 10 6 5 t D . (b) t R = 10 6 t D . (c) t R = 
10 5 r 't D . (d) t R = W 5 t D . (c) t R = 10 4 r 't D . (f) t R = 10 4 t D . 
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Figure 11. The L 2 distance Df(x,y,t) between f(x,y,t) and 
a near-invariant fi nv (x,y) computed for the second ensemble in 
FIG. 1, now perturbed by friction and white noise with = —E 
and variable r). (a) t R = W 6 - 5 t D . (b) t R = W 6 t D . (c) t R = 
10 5 S t D . (d) t R = 10 5 t D . (c) t R = 10 4 - 5 t D . (f) t R = 10 4 t D . 
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Figure 12. (a) A, the rate of approach towards a near-invariant 
distribution for the first ensemble of initial conditions in FIG. 1, 
allowing for friction and white noise with variable t R = -q^ 1 . (b) 
The same for the second ensemble, (c) T.05, the time required for 
the coarse-grained / associated with the first ensemble to evolve 
towards the near- invariant f„i v at the 5% level, (d) The same for 
the second ensemble. In each panel, the dashed line represents the 
value of A or r.os observed in the absence of friction and noise. 



a or t c , that the rate A associated with this exponential 
approach again scales logarithmically with r) or t R . 

As was observed also for other potentials (cf. Pogorelov 
& Kandrup 1999, Kandrup, Pogorelov & Sideris 2000), the 
response of an orbit ensemble seems insensitive to most de- 
tails. For fixed t c the two types of noise that were consid- 
ered had very similar effects; and the presence or absence of 
dynamical friction also proved largely irrelevant. All that re- 
ally seems to matter are the values of t c , which sets the time 
scale on which the noise changes appreciably, and t R , which 
probes the overall amplitude of the noise. In this sense, the 
results of these experiments were all consistent with the 
interpretation presented by Pogorelov & Kandrup (1999), 
namely that noise acts through a resonant coupling with an 
orbit. 

The spectral density, S(u>), given as the Fourier trans- 
form of the autocorrelation function K(t), characterises the 
degree to which the noise has significant power at different 
frequencies uj. The crucial point then is that the noise will 
have a significant effect on an orbit if and only if S has signif- 
icant power at frequencies co corresponding to the frequen- 
cies f2 ~ for which the orbit itself has significant power. 
It follows that, for autocorrelation times t c ^> to, coloured 
noise becomes comparatively ineffectual as a source of ac- 
celerated phase space transport. 

Examples of the effects of varying t c for fixed t R and 
fixed form of the coloured noise are illustrated in FIGURE 
13, which was generated once again from the ensembles of 
initial conditions used to generate FIGURE 8. The top two 
panels exhibit the effects of Ornstein-Uhlenbeck noise with 
tn, = 4000to on the convergence rate A; the middle two ex- 
hibit the effects of noise with autocorrelation function given 
by eq. (16). The bottom two panels exhibit r.os for the en- 
sembles used to generate panels (a) and (b). In each case, 
it is clear that, for t c <tc, the value of t c is essentially 
irrelevant and A assumes the value appropriate for white 
noise. Alternatively, for t c ^> to the noise has a compar- 
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Figure 13. (a) A, the rate of approach towards a near-invariant 
distribution for the first ensemble of initial conditions in FIG. 1, 
allowing for friction and Ornstcin-Uhlcnbeck coloured noise with 
tij = 4000trj> and variable t c . The dashed line corresponds to the 
rate A for unperturbed evolution, (b) The same for the second 
ensemble, (c) and (d) The same as (a) and (b), now allowing for 
coloured noise with autocorrelation function given by eq. (16). (e) 
t.o5, the time required for the coarse-grained / associated with 
the first ensemble to evolve in the presence of Ornstein-Uhlenbeck 
noise towards the near-invariant /„j„ at the 5% level, (f ) The same 
for the second ensemble. 

atively minimal effect. For values of t c comparable to or 
somewhat larger than to, A is a smoothly decreasing func- 
tion which exhibits a roughly logarithmic dependence on t c . 

That A only begins decreasing significantly for some- 
what larger values of t c for the fourth order noise with auto- 
correlation function given by eq. (16) reflects the fact that, 
for fixed tn, the autocorrelation function decreases some- 
what more slowly with t c than for the case of Ornstein- 
Uhlenbeck noise. This is, e.g. illustrated in FIGURE 14, 
which exhibits graphically the autocorrelation functions (12) 
and (15), in each case allowing for an autocorrelation time 
t c — 1.0 and a diffusion constant D = 1000. A(t c ) should 
not begin to decrease until t c becomes sufficiently large that 
there is an appreciable decrease in power for frequencies 
~ tp 1 relative to t c = 0. 

Interestingly, however, even very low frequency noise, 
with t c ^> to, can have an appreciable effect on quantities 
like T.ob- Even though such low frequency noise does not 
dramatically accelerate the initial exponential approach to- 
wards a near-invariant distribution, it does decrease signifi- 
cantly the time required for the second ensemble to approach 
a near- invariant distribution. This fact, illustrated in panels 
(e) and (f), reflects the fact that, even though the noise does 
not significantly impact the "average" behaviour of quanti- 
ties like the dispersion a x , it does suppress the large ampli- 
tude fluctations which are evident in the right hand panels 
of FIGURE 9. 




0246810 024 6 810 

t t 



Figure 14. (a) The autocorrrelation function K(t) for an 
Ornstcin-Uhlcnbeck process, plotted for t c = 1.0 and D = 1000. 
(b) The same for the autocorrelation function given by eq. (16). 

4 SUMMARY 

The aim of this paper, along with Paper I, has been to in- 
vestigate the role of chaotic processes in the triaxial Dehnen 
potential, which has been considered a realistic approxima- 
tion for at least the central regions of many elliptical galaxies 
and, possibly, also for some spiral galaxy bulges. Particular 
emphasis was placed on the way in which these processes 
are affected by the presence of internal and/or external per- 
turbations, modeled as periodic driving (Paper I) and/or 
dynamical friction plus white or coloured noise. Such per- 
turbations should in fact be operative in most real galaxies. 
The major findings can be summarised as follows: 

(i) In the presence of a central density cusp (7 > 0) 
and/or a central point singularity (black hole), the frac- 
tion / of 'strongly chaotic' orbits generally increases with 
increasing 7, with decreasing distance from the centre, with 
increasing deviation from sphericality towards axisymmetry, 
and with increasing deviation from axisymmetry towards tri- 
axiality. Typical values for / in triaxial configurations range 
between 30 and 70 percent. For a given deviation from spher- 
icality towards axisymmetry, / is higher when the deviation 
is prolate than when it is oblate. 

(ii) The fact that there are unquestionably regular orbits 
passing very close to the centre suggests that chaos is prob- 
ably triggered by appropriate resonance overlaps, which are 
stronger in steeper cusps, rather than 'close encounters' with 
the central cusp. 

(iii) In the absence of a central density cusp (7 = 0) and 
of a central black hole, the fraction of strongly chaotic orbits 
increases with distance from the centre. This is interpreted 
as evidence that the central cusp plays a dominant role in 
generating chaos for 7 > 0. 

(iv) The exponential instability of the strongly chaotic 
orbits, as measured by the value of their maximal Lyapunov 
exponent, is characterised by a timescale of order a few dy- 
namical times. However, unlike many other two- and three- 
dimensional potentials, chaotic orbits in the triaxial Dehnen 
potential tend to be extremely 'sticky,' even for integration 
times of 20, 000 to or longer. 

(v) The presence of a central black hole as large as 1% 
of the total mass causes a relatively small increase in the 
fraction of strongly chaotic orbits, mainly near the centre. 
However, the values of maximal Lyapunov exponents can 
increase substantially. 

(vi) Analogous results were also obtained for a toy model 
consisting of an anisotropic oscillator with a softened Plum- 
mer sphere (supermassive black hole) superimposed in the 
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centre (Kandrup & Sideris 2002). This suggests that they 
may be generic to nonaxisymmetric systems with central 
density cusps and/or central black holes. 

(vii) The evolution of an initially localised ensemble of 
chaotic orbits is typically characterised by three time scales, 
two comparatively short and the third much longer. The 
initial, shortest time scale evolution corresponds to an ex- 
ponential divergence at a rate A comparable to a typical 
value of the largest finite-time Lyapunov exponent \ f° r the 
orbits in the ensemble. This is followed by an exponential ap- 
proach towards a near- invariant distribution f n i v at a some- 
what slower, but still comparable, rate A Si A ~ \. This 
comparatively rapid evolution is then followed by a much 
slower evolution as the orbits diffuse along Arnold webs so 
as to, eventually, sample the true equilibrium. Since A -1 is 
typically comparable to a typical dynamical time, the first 
two stages correspond to timescales r ~ to', the time scale 
for the third stage is typically much longer, t ^> to- 

(viii) All three of these timescales can be shortened con- 
siderably by allowing for discreteness effects associated with 
two-body relaxation and modeled as a superposition of fric- 
tion and additive or multiplicative homogeneous Gaussian 
white noise in the context of a Langevin, or Fokker-Planck, 
description. 

(ix) Coloured noise, corresponding to finite-duration ran- 
dom kicks characterised by an autocorrelation time t c , can 
have the same effects as white noise, provided that t c Z^to- 



5 IMPLICATIONS FOR REAL GALAXIES 

In assessing the implications of these results for real galaxies, 
it is useful to differentiate between three relatively distinct 
scenaria, which are discussed in more detail in the following 
sections: 

(i) No external perturbations: Relatively isolated 'field' 
galaxies, for which the time-independent Hamiltonian de- 
scription is largely adequate. However, internal discreteness 
effects are still at work, and can affect the efficiency with 
which such a galaxy can reach a (quasi-)equilibrium. 

(ii) Weak external perturbations: Galaxies with a small 
number of satellite/companion objects and/or belonging to 
a moderately dense cluster of galaxies, where the (internal) 
galactic bulk potential still remains approximately time in- 
dependent. Internal perturbations can still be important, as 
in isolated galaxies. 

(iii) Strong perturbations: The later stages of violent in- 
teractions where the time-independent approximation com- 
pletely fails. Examples of such interactions include galaxy 
mergers and frequent high-speed galaxy encounters ('galaxy 
harassment'). In such cases, a 'violent relaxation' approxi- 
mation is more appropriate. 

This discussion does not include the effects of a dissipa- 
tive component (cold gas), so it is perhaps most appropriate 
for elliptical galaxies and for certain spiral bulges. 



5.1 Isolated galaxies (no external perturbations) 

Important discreteness effects inside a galaxy include: 



(i) Stellar (two-body) encounters, where the distance be- 
tween objects is very small compared with the size of the 
galaxy. These can be modeled adequately as instantaneous 
(t c <C to; cf. Eq. (17) in Section 5.2) random kicks, i.e., as 
friction and white noise characterised by a relaxation time 
t R = n' 1 and a temperature 9 = -E. FIGURES 10 and If 
exhibit the effects of white noise for a few values of E and 
several values of tR. 

(ii) Large-scale organised motions and encounters with 
massive, extended objects such as giant molecular clouds 
and star clusters. For these somewhat larger objects, t c is 
no longer completely negligible compared with to- However, 
to the extent that t c /tD ^ 0.2 or so, the results of Section 3.4 
indicate that white noise should still be adequate. Whether 
or not this is the case can depend on details such as the ve- 
locity field in the galaxy. For instance, t c /tD can be <C 1 in 
'thermally' supported objects, such as many elliptical galax- 
ies, where the duration of encounters can be short. However, 
in locales where the velocity dispersions are small and mo- 
tion is highly organised, such as in disks of spiral galaxies, 
it may be that t c ~ to or even t c 3> to (consider, e.g., the 
long 'encounter' of a star with a molecular cloud, both mov- 
ing along near-circular orbits of similar radii in the disk of a 
spiral galaxy). Such encounters can be seen as random kicks 
of finite duration which could be modeled as friction and 
coloured noise with appropriate autocorrelation times t c . 

Although the presence and properties of this coloured 
noise depend on the details of the internal kinematics of 
the galaxy, its dynamical consequences should be insensi- 
tive to these details. As shown in Section 3.4, the effects of 
coloured noise are much the same as for white noise when 
t c S2 to; and since, for fixed amplitude, coloured noise has 
a much weaker effect for t c 3> to, it should prove largely 
irrelevant for isolated galaxies which do not experience the 
effects of very massive extended objects. 

What are the dynamical implications of these inter- 
nal irregularities? Conventional wisdom (cf. e.g. Binney & 
Tremaine 1987) assumes that the dynamics of isolated galax- 
ies is governed solely by the bulk potential of the constituent 
stars because discreteness effects such as two-body interac- 
tions do not alter appreciably the values of the integrals of 
motion (such as the energy E) over the course of a Hubble 
time (cf. Chandrasekhar 1941). This treatment should, in- 
deed, be valid for integrable potentials where the motion is 
completely regular, as well as for hyperbolic nonintegrable 
potentials where motion is completely ergodic. 

However, it may not be appropriate to neglect discrete- 
ness effects in potentials admitting a large measure of both 
regular and chaotic orbits, such as the triaxial Dehnen po- 
tential and the majority of nonspherical potentials used in 
galactic dynamics, even though the values of the integrals of 
the motion still do remain approximately constant. Chaotic 
orbits in such potentials are often 'sticky,' i.e., when inte- 
grated over some finite time, they can spend a large frac- 
tion of that time 'trapped' inside cantori or within distinct 
phase space regions separated by the Arnold web. However, 
as explained in the preceding sections, discreteness effects 
can dramatically accelerate the rate at which sticky orbits 
become 'untrapped' and approach their respective (near-) 
invariant distributions. 
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These considerations can be important in at least two 
ways. First, discreteness effects, by acting as a continuous 
'microharassment' that violates Liouville's theorem, tend to 
smooth away substructures of phase space which correspond 
to regions of extreme instability, such as homoclinic points or 
high-order cantori. This can be of relevance when using fixed 
potentials as realistic galactic models. A study of the motion 
near unstable regions in such potentials should probably use 
realistic levels of noise during orbit integrations, correspond- 
ing to the graininess expected in the actual galactic poten- 
tial, in order to improve the robustness against structural 
perturbations in the form and in the Hamiltonian character 
of the potential. 

A second consideration involves the construction of 
self-consistent equilibrium models of galaxies using or- 
bital superposition methods such as Schwarzschild's method 
(Schwarzschild, 1979) and more recent variants thereof. The 
objective of such methods is to find a weighted superpo- 
sition of orbits that reproduces the mass distribution gen- 
crating the potential in which the orbits were evolved. In 
order for the model to be a true equilibrium, each su- 
perimposed ('library') orbit should correspond to a time- 
independent building block, i.e., to a time-averaged (and 
hence space-averaged) approximation of the invariant dis- 
tribution formed by galaxy orbits which share a particu- 
lar combination of values for the applicable integral(s) of 
the motion. However, a library orbit which remains sticky 
throughout much of its numerical integration time may not 
be a good time-independent building block since it does not 
sample uniformly its invariant measure. This remains true 
even if there are 'real' sticky orbits, in the 'real' galaxy, 
which never become 'unstuck' over a Hubble time! It is, 
therefore, desirable to integrate library orbits using noise 
with an amplitude chosen so that it accelerates phase space 
transport as much as possible, without adversely affect- 
ing the properties of the potential. In practice, this usually 
means that one can safely use a noise amplitude correspond- 
ing to the expected relaxation time of the galaxy or the ac- 
curacy of the numerical integrator, whichever is greater. If 
this level of noise is not enough to disentangle the sticky or- 
bits, then one could experiment with further increasing the 
amplitude, as long as the salient properties of the potential 
remain unaffected, and/or increasing the integration time. 

It should be noted at this point that it is not always re- 
alistic to expect that the invariant measure can be sampled 
well for all potentials using integration times which are not 
excessively long. The extreme stickiness exhibited, for ex- 
ample, by the cuspy triaxial potentials studied in this paper 
and in Paper I, makes it very hard to achieve a good nu- 
merical approximation of the invariant measure. In a strict 
sense, this would mean that it is not possible to construct 
equilibria of such potentials using Schwarzschild-like meth- 
ods. However, it may also be the case that 'partially mixed' 
near-equilibria can be realistic approximations over a Hub- 
ble time, especially when the partial mixing refers mostly 
to the outer parts of the galaxy, where tD is very long. One 
cannot properly answer this question without actually con- 
structing Schwarzschild models for these potentials. 

Finally, a legitimate question from a numerical stand- 
point concerns the properties of the noise that should be 
used (other than its amplitude), and whether one should 
use white or coloured noise or some combination of the two. 



The good news here, as explained in the preceding sections, 
is that the details of the noise are largely irrelevant, as long 
as the autocorrelation time t c to- Since white noise can 
be seen as coloured noise with t c — » 0, and considering that 
coloured noise is harder to code and requires considerably 
more CPU cycles to compute, it follows that simply us- 
ing white noise of the appropriate amplitude is completely 
adequate in most cases. It is interesting to note that the 
'numerical noise' associated with the integration of orbits 
could, perhaps, play the role of a white noise generator, un- 
der appropriate circumstances; however, this is a question 
that would require a separate investigation. 



5.2 Weak external perturbations 

The discussion in the preceding subsection concerning the 
consequences of discreteness effects in isolated galaxies, ob- 
viously applies fully also in the case of a galaxy surrounded 
by other objects. However, these objects act as additional 
sources of time-dependent perturbations, and their effects 
are inversigated here. 

In a first approximation, the effects of a companion 
and other nearby objects can be modeled in the spirit of 
Chandrasekhar's (1941) 'nearest neighbour approximation,' 
as justified by Chandrasekhar & von Neuman (1942) (see 
also Kandrup 1981), which associates the random gravita- 
tional force acting on any given object with one or two par- 
ticularly proximate neighbouring objects. Let v denote the 
typical speed of stars in the original galaxy and u the typi- 
cal speed with which other nearby galaxies are moving with 
respect to that galaxy. Similarly, let r be a measure of the 
physical size of the original galaxy and d the typical sepa- 
ration between galaxies. And finally, let M denote the mass 
of the original galaxy and m the mass of a typical nearby 
galaxy. 

The natural time scale on which the external forces 
acting on the original galaxy change significantly should 
be smaller than, or comparable to, the time required for 
a nearby galaxy to travel a distance ~ d, so that t c should 
be no larger than 



u \r \u , 



(17) 



where to ~ r/v denotes a characteristic dynamical time for 
the original galaxy. The force per unit mass associated with 
such a nearby galaxy will have a typical size 



F 



Gm 
d 2 



GM 
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(18) 



On dimensional grounds, tn satisfies (cf. Eq. (11)) 

(19) 



F 2 t c ~ v 2 t _1 



so that 



t R 



v 6 r 6 M\ d\ I u 



G 2 M 2 \ m 



(20) 



However, an application of the the Virial Theorem, 
GM/r ~ v 2 , then implies that 



(21) 
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Figure 15. (a) The L 2 distance Df(x, y, t) between f(x, y, t) and 
a near-invariant f n iv(^,y) for the first ensemble in FIG. 1, now 
perturbed by friction and Ornstein-Uhlcnbcck colored noise with 
© = — E, t R = 10 6 t£>, and t c = 4trj. (b) The same for the second 
ensemble in FIG. 1. (c) The same as (a) but with t R = W 5 t D . 
(d) The same as (b) but with t R = 10 5 t_o. (e) The same as (a) 
but with t R = V) 4 tu- (f) The same as (b) but with t R = 10 4 t_o. 
(g) The same as (a) but with t R = 10 3 to- (h) The same as (b) 
but with t R = 10 3 t D 



The presence of satellite/companion objects or of a 
nearby /surrounding cluster of galaxies will typically incur 
comparatively low-amplitude irregularities in the bulk po- 
tential associated, e.g., with a close encounter that has dis- 
placed the galaxy from a near-equilibrium. If such a system- 
atic time-dependence is present, t c could become apprecia- 
ble. In this case, the lowest frequency contributions might 
be expected to correspond to one or two quasinormal modes 
characterised by frequencies ui ~ t n 1 , but one might also ex- 
pect a larger collection of higher frequency modes which, in 
a first approximation, could be modeled as inducing random 
irregularities with characteristic time scale t c comparable to, 
but somewhat smaller than, to- 

Simple, but not unrealistic, models allowing for the 
effects of multiple satellite galaxies and other compan- 
ion objects might entail the choices d ~ 4r and u ~ v, for 
which t c ~ Ato and t R ~ 64(M /rafto- Similarly, allowing 
for interactions between galaxies of comparable size in a 
dense environment such as provided by the Coma Cluster 
might entail d ~ 12r and u ~ v, for which t c ~ 12to and 
t R ~ 1728(M/m) to- Internal irregularities associated with 
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Figure 16. (a) The L? distance Df(x, y, t) between f(x, y, t) and 
a near-invariant fniv(z,y) for the first ensemble in FIG. 1, now 
perturbed by friction and Ornstein-Uhlcnbcck colored noise with 
= —E, t R = 10 4 t£>, and t c = Yltjj. (b) The same for the second 
ensemble in FIG. 1. (c) The same as (a) but with t R = 10 3 t c (c) 
The same as (d) but with t R = 10 3 id 
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Figure 17. (a) The L 2 distance Df(x,y,t) between f(x,y,t) 
and a near-invariant f n iv(%,y) for the first ensemble in FIG. 1, 
now perturbed by friction and Ornstein-Uhlcnbcck colored noise 
with = -E, t R = 10 6 t D , and t c = 0At D . (b) The same for 
the second ensemble in FIG. 1. (c) The same as (a) but with 
t R = 10 5 t D . (d) The same as (b) but with t R = 10 s t D . 



changes in the bulk potential might entail t c ~ 0.4to and a 
large range of values for t R . 

FIGURE 15 exhibits the effects of Ornstein-Uhlenbeck 
noise with t c = 4t_o and variable t R ranging from 10 3 tu and 
10 6 . FIGURE 16 exhibits the effects of Ornstein-Uhlenbeck 
noise with t c = 12to and t R = W^to and 10 4 trj. It is clear 
that, for t c — 4to, t R — 10 5 to has an appreciable and 10 4 io 
a large effect. Similarly, for t c = 12to, a relaxation time 
t R — 10 4 tD is short enough to be important and 10 3 to again 
has a large effect. It thus follows that other giant galaxies 
with m ~ M may be expected to have significant effects in a 
dense cluster environment where d is as small as ~ lOr, and 
that companion/satellite galaxies with m ~ 0.1M or even 
smaller could also play an important role. 

FIGURE 17 exhibits the effects of Ornstein-Uhlenbeck 
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noise with t c = QAto and variable tR. Here it is clear that 
for tn as small as 10 6 to the noise can have an apprecia- 
ble effect. This would suggest that comparatively low level 
time-dependent irregularities associated with a galaxy out 
of equilibrium could indeed play a significant role in accel- 
erating the approach towards a near-equilibrium. 

What are the dynamical consequences of these exter- 
nal perturbations? Even though they are characterised by 
timescales which are long compared with to they can have 
a substantial effect. For fixed amplitude, coloured noise with 
autocorrelation time t c to will have a much stronger effect 
than noise with t c 3> to- However, external perturbations 
associated with nearby massive objects will in general have 
much larger amplitude and, at least in the outer portions 
of the galaxy, t c may not be that much longer than to- 
Perturbations with t c long compared with to can still be 
important as an evolutionary mechanism if the amplitude 
is sufficiently large. This would suggest that, even if viable 
for isolated galaxies, 'partially mixed' equilibria are not an 
option in high density environments! 

5.3 Strong perturbations (violent relaxation) 

The discussion hitherto has focused on orbit ensembles 
evolved in a fixed potential, possibly perturbed by low am- 
plitude perturbations. However, chaotic phase mixing also 
has important potential implications for the evolution of 
systems that evidence a strong time dependence, such as 
galaxies in the process of merging. Indeed, one might argue 
(e.g., Kandrup, Vass, and Sideris 2003) that chaotic phase 
mixing is an important component which must be incorpo- 
rated into any successful theory of violent relaxation. Ar- 
guably, the crucial ingredient of Lynden-BelPs (1967) origi- 
nal proposal is that phase mixing induces a coarse-grained 
approach towards equilibrium, his prototypical example be- 
ing the motion of particles in a one-dimensional 'pig-trough,' 
in which all the orbits are regular. The problem, however, 
with this and other similar examples is that the approach 
towards equilibrium is much less efficient than what simu- 
lations suggest for real stellar systems (cf. Kandrup 1999), 
typically proceeding at best as a power law in time. If, how- 
ever, one considers a flow that is chaotic rather than regular, 
the approach towards equilibrium should proceed exponen- 
tially in time, which implies that a near-equilibrium can be 
achieved much more quickly. 

One might perhaps object to this argument on the 
grounds that realistic systems will in general contain sig- 
nificant numbers of both regular and chaotic orbits, and 
that chaotic phase mixing should dominate the evolution 
of a time-dependent system only if chaotic orbits are much 
more common than regular orbits. Otherwise chaotic phase 
mixing would not be sufficiently ubiquitous as to have a 
dramatic effect on the system as a whole. It is certainly 
likely that, if only a small fraction of the orbits are chaotic, 
chaotic phase mixing will only have a comparatively minor 
effect. However, there is solid reason to believe that time- 
dependent potentials tend to admit much larger numbers of 
chaotic orbits, especially if the potential exhibits pulsations. 

Theoretically it is easy to understand why this might be 
so. If an orbit is regular, it must be restricted by constants of 
the motion, either global integrals like angular momentum 
associated with global symmetries, or 'local' integrals (cf. 



Lichtenberg & Lieberman 1992) which restrict the motion of 
regular orbits in nonintegrable systems. Making the poten- 
tial time-dependent removes the symmetry associated with 
time translation invariance so that neither the energy nor the 
Jacobi integral is conserved. In a generic time-independent 
potential, regularity requires two local integrals; in a generic 
time-dependent system three local integrals are required. 
One might thus anticipate that, at any given time, a larger 
fraction of the orbits will be chaotic, provided at least that 
the time-dependence is sufficiently strong that the energy 
cannot be treated as (nearly) an adiabatic invariant. 

For the case in which the time-dependence involves 
large amplitude systematic oscillations, one can in fact 
demonstrate that, in many cases - probably generically - the 
relative measure of chaotic orbits and the size of the largest 
finite time Lyapunov exponent will both increase dramati- 
cally (Kandrup, Vass, and Sideris 2003). Indeed, it has been 
recognised in both the accelerator dynamics (cf. Gluckstern 
1994) and the nonneutral plasma communities (cf. Stras- 
burg and Davidson 2000) that chaos associated with such a 
resonance can be important at a practical level. In partic- 
ular, it can result in particles in the core of an accelerator 
beam being ejected into an outerlying halo, thus resulting 
in a highly undesirable increase in the size of the beam. 

The basic idea is that the introduction of an oscilla- 
tory time-dependence can trigger a parametric resonance, 
involving a coupling between the frequencies associated with 
the pulsations and the frequencies of the orbits, making 
many/most of the orbits chaotic. The important point then 
is that, even for moderate amplitude oscillations, this reso- 
nance can be very broad, requiring only that the two sets 
of frequencies be comparable to within an order of magni- 
tude or so. Given, however, that there is only one natural 
time scale for the problem, namely the dynamical to, which 
sets both the characteristic orbital time scale and the pulsa- 
tion time, one would anticipate that, throughout most of the 
galaxy, the frequencies will be sufficiently close to trigger the 
resonance. Simple toy models involving an integrable Plum- 
mer sphere subjected to damped oscillations can, within 
a time ~ 10to, exhibit both near-complete chaotic phase 
mixing and an approach towards a nearly time-independent 
state (Kandrup, Vass, and Sideris 2003). 

The same physical processes can also act on shorter 
scales. A supermassive black hole binary introduced into 
the center of a galaxy - with or without a cusp - serves 
as a time-dependent perturbation which can trigger large 
amounts of chaos even at distances from the center much 
larger than the binary orbit; and the resulting chaos can lead 
to efficient phase mixing which occasions significant changes 
in the density distribution and, hence, the observable sur- 
face brightness profile (Kandrup, Sideris, Terzic, and Bohn 
2003). 

The work described here and in Paper I reinforces the 
expectation that chaos may be ubiquituous in cuspy galax- 
ies, especially those that manifest significant deviations from 
axisymmetry, and that that chaos may have significant im- 
plications for both structure and evolution. In particular, 're- 
alistic' galaxies, characterised by a bulk potential that admit 
a complex coexistence of regular and chaotic behaviour, may 
be substantially more susceptible to various 'irregularities' 
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than had been recognised originally. The crucial remaining 
issue, which has yet to be addressed, would seem to be pre- 
cisely how the effects described here will manifested in the 
context of a time-dependent, fully self-consistent evolution. 

Finally, however, one can conclude by comparing the re- 
sults derived in this paper with models which might naively 
seem to contradict one of its principal conclusions (as well as 
of Gerhard & Binney [1985], Merritt & Fridman [1996], Mer- 
ritt & Valluri [1996], and Merritt [1997]), namely that chaos 
should be ubiquituous in cuspy triaxial galaxies. Specifically, 
Holley-Bockelmann et al (2001) (hereafter HB) constructed 
JV-body realisations of triaxial (quasi-) equilibria with a cen- 
tral density cusp by adiabatically 'squeezing' spherical 7 = 1 
Hernquist models and, for one particular choice of moder- 
ately triaxial axis ratios (a : b : c = 1 : 0.85 : 0.7) found no 
chaotic, or at least no strongly chaotic, orbits. And similarly, 
Poon & Merritt (2002) (hereafter PM) were able to construct 
Schwarzschild models with 7 = 1 and 2 profiles containing 
only regular orbits (although they also constructed models 
with both regular and chaotic orbits). iV-body realisations 
of all their Schwarzschild models revealed only small fluctu- 
ations over a few crossing times. 

In point of fact it is not completely clear that there 
are no chaotic orbits in the HB model. As those authors 
themselves stated, their algorithm to identify chaotic orbits, 
based on their Fourier transforms, could easily have missed 
a large number of nearly-regular, but still chaotic, orbits. 
Moreover, since their search for chaos involved integrating 
initial conditions in the fixed iV-body background associated 
with a single t = constant snapshot, they excluded explicitly 
the possibility that weak amplitude oscillations which, as 
is evident from Figure 3 of HB, exist in the real model, 
could trigger chaos via parametric resonance. This model 
is, however, clearly interesting in that it may represent an 
example of a near-equilibrium supported using only 'local' 
(i.e., 'third') integrals. All three of these points can also be 
made for the PM models. 

In any event, it could well be that neither the triaxial 
Dehnen models nor the HB or PM models constitute true 
time- independent self-consistent equilibria. However, as em- 
phasised in Paper I (and discussed more extensively in Kan- 
drup [2002]), it is quite possible that a galaxy will settle 
initially into one triaxial near-equilibrium state and then 
evolve slowly through a sequence of near-equilibria with- 
out necessarily becoming more axisymmetric. There is at 
the present time no compelling evidence that a slow secular 
evolution necessarily involves an evolution towards axisym- 
metry. What might, however, be true is that a galaxy that 
resembles a triaxial Dehnen model with large measures of 
wildly chaotic orbits would evolve towards less chaotic con- 
figurations better represented by an HB- or PM-type model. 
The work described in the present paper could provide po- 
tentially significant insights as to precisely how such an evo- 
lution might proceed. 
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